Collective Oscillations of Vortex Lattices in Rotating Bose-Einstein Condensates 
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The complete low-energy collective-excitation spectrum of vortex lattices is discussed for rotating 
Bose-Einstein condensates (BEC) by solving the Bogoliubov-de Gennes (BdG) equation, yielding, 
e.g., the Tkachenko mode recently observed at JILA. The totally symmetric subset of these modes 
includes the transverse shear, common longitudinal, and differential longitudinal modes. We also 
solve the time-dependent Gross-Pitaevskii (TDGP) equation to simulate the actual JILA experi- 
ment, obtaining the Tkachenko mode and identifying a pair of breathing modes. Combining both 
the BdG and TDGP approaches allows one to unambiguously identify every observed mode. 
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Owing to their fundamental significance for superflu- 
idity, quantized vortices have attracted widespread inter- 
est in many different physical systems ranging from su- 
perconductors, superfluid 3 He and 4 He, to neutron stars 
[l| and extending to cosmology. Recently, vortices have 
been created in dilute alkali-atom gases using three dif- 
ferent methods: phase imprinting 2], topological phase 
engineering Q and optical spoon stirring EL The for- 
mer two were first predicted theoretically [fj |g ■ Several 
groups are now able to routinely prepare a vortex array 
with hundreds of vortices in a BEC. 

In the mid 1960's, Tkachenko [7j predicted that a vor- 
tex lattice would sustain a collective vortex-oscillation 
mode in which the vortex cores move elliptically around 
the equilibrium positions. Subsequent theoretical ad- 
vances were made Q, 0, H, 0, ^J. In superfluid 4 He, 
the Tkachenko modes were first observed in 1982 fl2|| . 

Recently, Coddington et al. 01 succeeded in observ- 
ing the Tkachenko mode (TK) in a rotating BEC of 87 Rb. 
The TK wave was excited by removing condensate from 
the central region of the rotating cloud. They created 
the lowest and second-lowest TK modes and measured 
their energies u>\$ and W2,o as functions of the rotation 
frequency f2. They also discovered various new phenom- 
ena, some of which were explained by two groups [lfil,llfij |: 
Anglin and Crescimanno extended the previous hy- 
drodynamic description for infinite systems to a finite 
harmonically trapped system, and Baym discovered 
subtle effects due to the finite compressibility. 

Here we choose a different approach to this problem: 
We wish to treat the whole low-lying collective-excitation 
spectrum of trapped BECs; not only the TK mode but 
also the other important modes and their intrinsic rela- 
tionships. We construct a first-principles theory, namely 
the Bogoliubov-dc Gennes equation (BdG) coupled with 
the Gross-Pitaevskii equation (CP). The set of equations 
within the Bogoliubov framework is regarded as the fully 
microscopic theory of dilute Bose gases, in the sense that 



there remain no adjustable parameters once we have fixed 
the atomic species and the atomic number. Thus it is 
quite reasonable to expect that this formalism must well 
apply to analyzing the TK mode and also to provide 
the complete spectral features of the low-energy excita- 
tions, beyond any limitations of hydrodynamics which 
has thus far been the only way to describe the TK mode. 
We can then better characterize the various modes, such 
as the three classes of compressional modes: the trans- 
verse, common longitudinal and the differential longitudi- 
nal waves. These are characteristic to the two-component 
system consisting of a vortex lattice and the superfluid. 
Some of the selected common longitudinal modes in a 
vortex lattice (the breathing and quadrupole modes) have 
been recently examined 17] within hydrodynamics. This 
theory may help one to gain improved understanding of 
related problems, such as vortex pinning and vortex melt- 
ing in superconductors which have thus far only been 
analyzed phenomenologically. 

In a frame rotating with the frequency Q, the time- 
dependent Gross-Pitaevskii equation (TDGP) may be ex- 
pressed for the condensate wavefunction tp [18j as 
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where p(O) = —ihW — mtlzxr and the effective confining 
potential is \4ff(r, fi) = |m(^ — il 2 )r 2 with the radial 
trap frequency u> r . 

In order to study the collective oscillations of vortex 
lattices microscopically, we first consider the equation of 
motion for a small perturbation around the stationary 



state (f> g , i.e., ip(r, t) = <p g (r) + u q (r)e 



B!(r)e w i f , 



where the equilibrium state 4> g is determined by the sta- 
tionary GP equation. By retaining terms up to first order 
in u and v, we derive the BdG equation: 
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where £(r, fi) = p 2 (ft)/2m + y cff (r, O) - /i + 2 5 |0 s (r)| 2 . 
We consider the JILA experiment with 2.0 x 10 6 atoms 
of 87 Rb confined in a trap with the radial frequency 
u> r /2rt = 8.3 Hz and an axial one uj z /2it = 5.2 Hz |l3| . 
Employing Thomas-Fermi (TF) theory and the assump- 
tion of solid-body rotation, the condensate aspect ratio 
is given as Atf = -Rtf/^tf oc (lu 2 — fJ 2 ) -1 / 2 where Rtf 
and Ztf are the condensate lengths along the r- and 
z-axes. This relation allows one to assume the system 
to constitute a two-dimensional geometry at high rota- 
tion frequencies. Under this assumption, we introduce 
the linear density n z (SY) along the z axis. Therefore, the 
equilibrium state 4> g must fulfill the normalization con- 
dition n z (Q) — J J \4> g \ 2 dxdy, where the linear density is 
obtained as n z (Q) — R^ F /16adf, with a an s-wave scat- 
tering length, and d 2 = h/(m^uj 2 — f2 2 ). We discretize 
the two-dimensional space typically into a 300 2 ~ 1000 2 
mesh to solve the TDGP and BdG equations. 

Here, since we consider a vortex array with sixfold sym- 
metry, the wavefunction of the stationary state obeys the 
condition, <p g (R n r) = 4> g (r)e" l7T / 3 , where R n r describes a 
rotation nir/3 (n integer) around the center of the trap, 
i.e., R n r — (xcos (mr/3) — 2/sin(n7r/3),a;sin(n7r/3) + 
y cos (rnr/ 3)). We then obtain the following relation 
from Eq. J2J): w q . m (i? n r) 
and 



j(r)exp [i2E( m +l) 



' q ,m(-R n r) = w q (r) exp [^(m — 1)1 , where ,„ — 
0, ±1,±2,3. In order to classify the collective excita- 
tions, we introduce the following classifying function: 

F<?\m) = /dr Wq (r)[ELo^™(^)]//dr| Uq (r)| 2 , 

which tends to 1 for a suitable m and for the others. 
Furthermore, we define the average angular momentum 
as qg = [(L z ) u + (L z ) v - (L z )^]/ J dr[\u(r)\ 2 + \v(r)\% 



where (L z )^ 



J dr(j)*L z (f) g / J dr\<j> g \ 2 , and (L z 



In an axisymmetric situation, m and 
qg merge into the same integer quantum number. 

The low-energy excitations are illustrated in Fig. 1. 
Here, the equilibrium states are found numerically via 
imaginary time propagation of Eq. Q; t — * r = —it. 
The initial vortex configuration is taken as a regular ar- 
ray with sixfold symmetry. The resulting configuration at 
il = 0.7uj r is formed by 37 vortices, displayed in Fig. 1(a). 
In Fig. 1(b), the excitation energies up to w = 3.5w r 
are shown as functions of the average angular momen- 
tum qg. Each collective mode is classified by a symmetry 
index, ro, obtained from the above function, (m) , 
which characterizes the oscillation pattern. For example, 
the breathing (BR), dipole (DP), and quadrupole (QP) 
modes have m = 0, ±1, ±2, respectively. The branch ex- 
tending from the origin towards higher qg-values consists 
of surface modes which are spaced periodically with m 
(mod 6), where the outer condensate surface oscillates, 
the vortex cores being almost stationary at their equi- 
librium positions. The other branch situated at higher 
energies and starting at u> = 2uj r represents another class 
of surface modes characterized by a node along the radial 
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FIG. 1: (a) Density profile of the condensate with 37 vortices 
at O = 0.7cu r . (b) Collective excitations up to w = 3.5uv. 
(c) Lowest excitations marked by a dotted line as functions 
of qg are displayed as functions of the symmetry index m. 
Here, TK, BR, DP, and QP denote the Tkachenko, breathing, 
dipole, and quadrupole modes, respectively. 



direction. The energy spectra for the low-lying modes are 
depicted in Fig. 1(c) as functions of the symmetry index 
m. These low-lying eigenstates are the lattice-oscillation 
modes coupled with the condensate motion and lead to 
the distortion of both the lattice and the condensate sur- 
face towards an m-fold symmetric shape. In addition to 
the lowest and second-lowest TK with m = 0, embedded 
among the other modes, we found a parallel-precession 
mode with m = — 1 and having the lowest energy where 
all the vortices precess in phase. For increasing energy, 
the modes for each m feature nodes along r and/or 8. 
The TK mode (with m — 0) is selectively excited by the 
Gaussian laser beam with 0-fold symmetry in the experi- 
ment [l3T |. Likewise, a mode with m = ±£ may be excited 
bv a disturbance (e.g., magnetic) with £-fold symmetry 
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FIG. 2: Oscillation patterns for the (a) lowest TK, (b) second 
TK, (c) lowest BR, and (d) second BR modes are indicated 
by the filled circles and solid lines. Empty circles and dotted 
lines correspond to the equilibrium state, <fi g . 

Figure 2 describes the oscillation patterns for se- 
lected collective modes obtained from a linear fluctuation 
<5V>(r, t): The lowest TK mode is shown in Fig. 2(a) 
where the empty (filled) circles indicate the equilibrium 
(quarter-period) positions of the vortex cores. The line 
is a fitted sinusoidal modulation wave (for clarity, the 
amplitude is magnified by a factor of 30). It is ob- 
served that the modulation displays a node at r no d e — 
0.64i?TF, which is in fair agreement with the JILA data 
of r no dc = 0.665i?TF- The second TK mode W2,o is shown 
in Fig. 2(b). It is seen that in contrast to the first TK 
mode whose core modulation is almost purely transver- 
sal, the second TK contains a longitudinal component, in 
addition to the transverse one. Thus, the core executes 
elliptic motion around its equilibrium position. This be- 
havior coincides with the general trend for the elliptic- 
ity of a TK mode, namely the expression for the ampli- 
tude ratio of the longitudinal and transverse components 
<5l/<5t oc l/yT2fcTK, where /ctk is the wavelength of the 
TK mode 1]. Therefore, the second TK with a shorter 
/ctk exhibits more pronounced longitudinal motion. We 
also point out that via encountering the condensate ra- 
dius in Figs. 2(a) and 2(b), the TK modes tend to accom- 
pany the condensate motion. We also observe that the ra- 
tio of the energies for these TK modes is u>2,o/^>x,o = 1-56 
at Q = 0.7uj r , which compares favorably with the exper- 
imental data of 1.8 ± 0.2 at fl = 0.95w r [l3| and with 
the value 1.63 obtained by Anglin and Crescimanno [lfij . 
who assert that it would be independent of S7. 

As typical vortex-lattice oscillations apart from the 
TK, we now consider two contrasting compressional 
modes which belong to the totally symmetric (m = 0) 



"breathing" branch. Figures 2(c) and 2(d) depict the 
oscillation patterns for the lowest (1BR) and the second- 
lowest (2BR) breathing modes in Fig. 1(b). The for- 
mer (latter) BR mode executes mutual in-phase (out- 
of-phase) motion between the cores and the condensate 
surface, as indicated with the two parallel (opposite) ar- 
rows in the figure. The inner (outer) arrow denotes the 
motion of the cores (condensates). Therefore, the for- 
mer belongs to the common longitudinal mode and the 
latter to the differential longitudinal mode. It should be 
pointed out that the inner core motion in Figs. 2(c) and 
2(d) is predominantly transversal; consequently, even in 
the nominally longitudinal BR mode all the vortex mo- 
tion is not necessarily purely longitudinal. This may be 
related to the 's-bending' [13| of the third observed mode. 

Having obtained the complete spectrum of the low- 
lying excitations, let us now proceed to extend the analy- 
sis into nonlinear dynamics with real-time evolution sim- 
ulated with the TDGP Eq. (JTJ. Since Coddington et al. 
excite the TK mode by shining resonant laser light in 
the center of the BEC to produce an inward flow, we in- 
troduce a Gaussian potential V(r) = —Tioj r e~ r /( r tf/2) 
localized in the center of the trap for a certain period and 
follow the time evolution according to Eq. (1). 

We illustrate the resulting oscillation patterns and 
their Fourier analyses in Fig. 3. At £1 = 0.7w r , the 37 vor- 
tices shown in Fig. 1(a) form a concentric regular lattice 
around the central vortex, consisting of three circles, j = 
1,2,3. We analyze the transverse vortex motion Yj^ n (t) = 
{xj,n{t),yj iTl (t)) in terms of the averaged angle 0j(t) = 
\ E«=o [ aTcta ^(yj,n(t)/xj tn (t))-axctan(y jtn (0)/xj tTl (0))] 
where n denotes the vortices aligned along the line and 
extending the angle nir/3 from the vertical, see Fig. 1(a). 
It is apparent in Fig. 3(a) where we plot the time depen- 
dence of the vortex motion for each circle, j = 1,2,3, 
that (i) there exist several superposed oscillations, (ii) 
The outermost vortices (j — 3) are in opposite phase 
with the inner vortices (j — 1,2). This result coincides 
with the above calculations based on the BdG equations. 
In fact, the nodal positions obtained from both calcula- 
tions agree quite well. 

In Fig. 3(b), the Fourier analyses of these transverse 
oscillations, 0j(t), and the longitudinal motions, rj(t) — 
\ En=o l r j,n(i)|, are depicted. It is seen that the sharp 
peak at w = 0.16a;,- precisely coincides with the first TK 
mode ojifi identified above, see Fig. 1(c). The second 
peak at u> = 2.0oj r corresponds to the 1BR mode, be- 
longing to the common longitudinal modes. The third 
peak at u) ~ 3.1cj r may be identified as the 2BR mode, 
see Figs. 1(b) and 2(d), and as a differential longitudinal 
mode. These three collective modes closely match the 
observed characteristics. In particular, the third mode 
which was tentatively assigned by Cozzini et al. and Choi 
et al. independently, as a higher-order hydrody- 

namic mode, is now identified above. It should be noted 
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FIG. 4: Frequencies of three modes obtained from the BdG 
and TDGP approaches are compared with the JILA data 
(filled squares, diamonds, and triangles) as functions of Q. 
In the main panel the TK (open squares), in the inset both 
the first BR (open diamonds) and the second BR (open trian- 
gles). The solid line represents the dispersion relation in Ref. 



FIG. 3: (a) Transverse oscillations of vortices on three con- 
centric circles, j = 1,2,3, from the trap center, (b) Fourier 
analyses of these oscillation patterns, 0j(t), and the longitu- 
dinal vortex motion, r, (i). 

that the resonances of these three modes for a Gaussian 
potential have been numerically reproduced over a wide 
range of rotation rates fl = 0.7 ~ 0.92, corresponding to 
37 ~ 121 vortices. 

In Fig. 4, we plot the first TK energy u>i t o as a func- 
tion of O/uv together with the experimental data and 
a hydrodynamic prediction by Anglin and Crescimanno 
(lif. There prevails close quantitative overall agreement 
between our results and the experimental data. We em- 
phasize that our calculations contain no adjustable pa- 
rameters and also that our computations within the BdG 
and TDGP approaches agree within numerical accuracy 
~ 10 _3 w r below < 0.8w r . Calculations for larger rota- 
tion rates, where BdG cannot be feasible from a numeri- 
cal point of view, are done with TDGP, which enables us 
to extrapolate the BdG results to larger rotation rates. 
The inset in Fig. 4 describes the O dependence of the 
1BR and 2BR modes. It is known that the former is 
consistent with an earlier prediction by Pitaevskii and 
Rosch |2(j , who point out that the 2D BR mode features 
the universal eigenfrequency ui = 2.0a;,-. Our result re- 
produces this result and further explains the observation 
mentioned above that the second peak in Fig. 3(b) is in- 
deed the 1BR mode. As for the third peak in Fig. 3(b), 
previously identified as the 2BR mode, it is reckoned from 
this inset that the observed value 18.5 ± 0.3Hz (~ 2.2w r ) 
appears accountable, judging from the overall Q depen- 
dence towards f2 c - 

In summary, we have discussed the complete low- 
energy excitation spectrum in a vortex lattice by solving 



the BdG equations. The m = subset of these solutions 
includes the transverse shear, common longitudinal, and 
differential longitudinal modes. We have also succeeded 
in simulating the actual experimental results, identifying 
a new pair of modes. 
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During the preparation of this manuscript, we learned 
about two closely related preprints |2l| and ji^. The 
former (latter) presents BdG (TDGP) treatments of the 
TK mode. 
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